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Electron-phonon coupling plays a central role for time-dependent phenomena in condensed matter, 
for example in photo-excitation experiments. We use the continuous-time quantum Monte Carlo 
method to study the real-time evolution of charge and spin correlation functions of a Peierls insulator 
after a quench to a noninteracting Hamiltonian. This approach gives exact results, and fully takes 
into account quantum phonon effects without relying on a Hilbert space truncation. It is also free 
from a dynamical sign problem. Since thermalization is absent for a quench to an integrable model, 
we observe dephasing over a time scale set by the electronic hopping integral. Our exact results 
provide a benchmark for more realistic calculations, and may be directly applicable to experiments 
with cold atoms or trapped ions. 
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I. INTRODUCTION 

The enormous current interest in the time evolution or 
nonequilibrium behavior of quantum systems is driven 
by significant developments in the field of experimental 
physics. Of particular importance are time-resolved stud- 
ies of cold atoms in optical lattices following the seminal 
work by Greiner et a/.Pand the progress with femtosec- 
ond spectroscopy)^] including time-resolved photoemis- 
sion spectroscopy!^ These methods enable researchers to 
carry out pump-probe experiments and to study photo- 
excitation and photo-induced phase transitions!^"^ 

Theoretical studies are particularly difficult for quan- 
tum systems with strong correlations, as relevant for cold 
atoms and several classes of materials. Whereas much 
progress has been made in developing methods suitable 
to study correlated systems in equilibrium, the descrip- 
tion of time-dependent or nonequilibrium problems is ex- 
tremely challenging. As a result, it is currently not pos- 
sible to reliably describe experiments. Promising direc- 
tions of research include time-dependent density matrix 
renormalization group methods,^ extensions of qua ntum 
cluster methods such as dynamical mean-field theorjWiLL^ 
and the variational cluster approach,^ as well as exten- 
sions of Luttinger liquid theory.^ A review of theoretical 
methods for correlated systems has been given by Eck- 
stein et alW 

Whereas these ideas have successfully been applied to 
models of interacting spins, fermions or bosons, much less 
is known about phenomena related to electron-phonon 
interaction. The latter plays a crucial role since phonons 
are often the fastest channel available for thermalization 
via scattering.^ However, the infinite bosonic Hilbert 
space as well as the resulting retarded electron-electron 
interaction pose a particular challenge to theoretical 
methods. Whereas exact results have been obtained for 
a single polaronj22HHl a single bipolaronpSl and a single- 
molecule devicep^ lattice models with finite band fillings 
are much more challenging. Recent work along these lines 
makes use of the density matrix renormalization group 
on small systems,^ the so-called double-phonon cloud 



method pH] restricted phonon Hilbert spaces p2l as well as 
classical phonons 

A particularly interesting class of experiments involves 
photo-induced p hase tran sition across a metal-insulator 
Peierls transitionPEHIUM] A laser pulse excites the sys- 
tem out of the insulating, charge-ordered state, as visi- 
ble for example from the time dependence of the gap in 
the band structure (probed by time and angle-resolved 
photoemission spectroscopy or from the reflectivity^ 
After the photo-induced melting of charge order, the sys- 
tem relaxes back towards its original state. Such exper- 
iments also reveal the typically very different timescales 
for electron and ion dynamics. So far, a photo-induced 
phase transition of a Peierls system has been considered 
theoretically in Ref. 29 using a classical phonon approx- 
imation. Since nonlocal correlations are a defining fea- 
ture of charge ordered states, local approximations such 
as dynamical-mean field theory are of limited use. 

Motivated by the ongoing effort to improve the avail- 
able theoretical tools for the study of electron-phonon 
systems, we consider here the quantum quench from a 
Peierls insulator to noninteracting fermions by switching 
off the electron-phonon interaction. The model we use 
is the one-dimensional Holstein model. The time evo- 
lution of charge and spin correlation functions after the 
quench is studied with an extension of the continuous- 
time quantum Monte Carlo method.^! Although such a 
quench does not reflect the complex experimental situ- 
ations encountered in pump-probe experiments, it pro- 
vides us with a rare opportunity to obtain exact re- 
sults for the time evolution of an electron-phonon system. 
Our findings can serve as a nontrivial benchmark for fu- 
ture numerical calculations (for example with methods 
that rely on a truncation of the phonon Hilbert space), 
and may even be relevant for experimental realiz ations 
of ele ctron -phonon models in optical lattice: 



or ion 



traps 
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The paper is organized as follows. In Sec. ITT] we dis- 



cuss the model, in Sec. Ill we provide details about the 
method used, in Sec. |IV| we present our results, and 
Sec. fVl contains our conclusions. 
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II. MODEL 

The Peierls metal-insulator transition in one dimen- 
sion can be investigated in the framework of Holstein's 
molecular crystal model, 3 ? defined by the Hamiltonian 

ka i \ / 

-g(t)Y,Qi i^i-i) . 

i 

The first term, referred to as Hq in the following, de- 
scribes nearest-neighbor hopping of electrons with the 
usual tight-binding band structure e(fc) = — 2Jcos(fc); 
we have set the lattice constant equal to one, and take 
the hopping integral as the energy unit. The second and 
third terms correspond to the lattice degrees of freedom 
and the electron-phonon coupling. The lattice is de- 
scribed by means of harmonic oscillators with frequency 
cjo = \J K/M, displacement Qi and momentum Pj, and 
the coupling is of the density-displacement type proposed 
by HolsteinJ^ g(t) is the coupling constant. We use 
the usual definitions h.- La — c\ a c ia , hi — X^"-^, and 
11 = {hi). A useful dimensionless coupling parameter is 
given by the ratio A = 2e p /W , where e p = g 2 /2K is the 
polaron binding energy, and W = 4 J is the bare band- 
width. We consider a half filled band, i.e., (hi) = 1. 

As discussed in Sec. [TJ we study a quantum quench to a 
noninteracting Hamiltonian. This is achieved by a time- 
dependent coupling constant g(t). Explicitly, we take 

9 (i) = { 3 ' i = °' (2) 
yw [0, t > 0. w 

Because there is no coupling between electrons and the 
lattice after such a quench, the time evolution will be 
determined by the hopping term Hq only. 



III. METHOD 

The continuous-time quantum Monte Carlo (CTQMC) 
method is based on a weak-coupling Dyson expansion of 
the path integral for the partition function!^! For finite 
systems, all possible vertex contributions can be summed 
stochastically, thereby giving exact results with only sta- 
tistical errors. The CTQMC method has been used quite 
extensively to solve impurity problems in the context of 
dynamical mean- field theory, see Ref. 00] for a review, as 
well as for lattice simulations. 

Because the method is formulated in terms of an ac- 
tion rather than a Hamiltonian, it permits us to simu- 
late models with retarded interactions. Such interactions 
naturally arise as a result of electron-phonon interaction, 
and models such as Eq. ([I]) can be simulated by integrat- 
ing out the lattice degrees of freedom analytically.^ This 



approach h as be en applied to study spinless^ and spin- 
ful Holsteiri- 43 * 44 ! and Holstein-Hubbard models,^ as well 
as models with a nonlocal electron-phonon interaction.^ 
Numerical details can be found in Ref. 251 Because the 
method relies on an analytical integration over phonon 
degrees of freedom, lattice properties such as local dis- 
placements cannot be directly measured. 

Quantum Monte Carlo approaches to time-dependent 
problems are inherently limited by the dynamical sign 
problem related to the time evolution operator exp(iiJt). 
Nevertheless, progress has been made in applying such 
methods to simplified situations where the sign problem 
can be avoided or limited, including quantum quenches 
to noninteracting systemsp^and studies of (effective) im- 
purity problems. For the CTQMC method used here, 
the first application to a nonequilibrium problem in- 
volving phonons was achieved in Ref. 1481 and a review 
of previous work is found in Ref. |40] Recently, the 
weak-coupling CTQMC method was formulated on the 
Keldysh contour.^ 

In principle, the method of Ref. |47] can be extended 
to include a retarded electron-electron interaction on the 
Keldysh contour. The resulting algorithm would permit 
us to study rather general problems, but with the usual 
limitation to short times originating from the dynamical 
sign problemP^ Instead, we follow a simpler route suffi- 
cient to study quenches to noninteracting systems, based 
on the idea used by Goth and Assaad in the context of the 
auxiliary-field quantum Monte Carlo methodP^ If we are 
only interested in the time evolution after the electron- 
phonon coupling g has been switched off [see Eq. pi], in- 
teraction vertices are restricted to imaginary times on the 
Keldysh contour. The time dependence of observables is 
then fully determined by the noninteracting Hamiltonian 
H a , the first term in Eq. ([l]). 

The numerical method proceeds as follows. The inter- 
acting Holstein model ([T]) is simulated with the CTQMC 
method using standard Monte Carlo updates (addition 
and removal of vertices, flipping of Ising spins ) to ob- 
tain configurations reflecting the initial, correlated state 
at time t = 0. The interacting single-particle Green's 
function matrix at equal imaginary times, 

[G„(t = 0)] y = <c iCT ct CT ) , (3) 

is calculated fo r a g iven vertex configuration from the 
Dyson equation! 32 ! 41 ! From G a (t = 0), all relevant corre- 
lation functions and o bserv ables at t = can be obtained 
using Wick's theorem! 32 ! 41 ! For quenches of the form ([2]), 
expectation values at times t > follow from the time- 
propagated Green function 

G a (t) = e~' l6ot G a (t = 0) e ikat , (4) 

in combination with Wick's theorem. The additional nu- 
merical effort compared to a time-independent simulation 
amounts to the evaluation of Eq. Q and of observables 
of interest for each value of t. Results for each t are aver- 
aged over different Monte Carlo configurations. Whereas 



3 



equal-time observables can easily be calculated in this 
way, time-displaced correlation functions require a for- 
mulation on the Keldsyh contour. 



IV. RESULTS 

For a half-filled band, the Holstein model describes 
a transition from a metal to a Peierls insulator driven 
by the electron-phonon interaction. At T = 0, the 
Peierls state has long-rage charge order with wavevec- 
tor q = 2kp, and a gap in the single-particle spectrum! 4 ^! 
In addition, a spin gap presumably exists for any nonzero 
electron-phonon interaction, as a result of electrons pair- 
ing into singlets. 45 Whereas the model is dimerized and 
insulating for an y A > in the mean-field limit wq = 
numerical resulta^HSU indicate a nonzero critical value for 
cjo > as a result of quantum lattice fluctuations. For 
an overview of previous work, see Ref. H5] 

Here, we consider parameters for which Hamilto- 
nian (!]) is in the insulating Peierls phase. Explicitly, 
we take A = 0.5, and different values of the phonon 
frequency, ujq = 0.5J, and 5J. According to quantum 
Monte Carlo resultspl A = 0.5 is well inside the ordered 
phase for uig = 0.5J (the critical coupling is about 0.25), 
and just above the critical coupling for ujq — 5 J. Since 
our simulations are carried out at a low but finite tem- 
perature (we consider a fixed system size L = 22 and 
pj = 2L = 44), the system will be in a thermal ini- 
tial state before the quench. However, the results shown 
below are practically identical for j3J = L and f3J = 2L. 

Among the equal-time observables accessible with the 
present method, the real-space charge and spin correla- 
tion functions 



S p (r) = (n r h ) , 
S a (r) — (S r S n ) , 



(5) 



are of particular interest. Since we consider a quench to 
the noninteracting Hamiltonian, we compare the corre- 
lation functions to those of free fermions. For the latter, 
Luttinger liquid theory gives^ 



S p (r) = S a (r) = - 



1 



C 



(7rr) 2 



cos(2fcFr) . 



(6) 



Here, we have set K p — K„ = I 



as appropriate for 
free fermions. In contrast, the Peierls state is formally 
described by K p — (implying long-range 2fcp charge 
order) and K a = (reflecting the spin gap which leads 
to exponentially suppressed spin correlations).^ 

Figure [l] shows results for S p (r,t) for ll>o — 0.5 J and 
A = 0.5, as calculated from the time-dependent Green's 
function Q. Panels (a)-(d) correspond to different times 
t in the interval [0,0. 7 J], whereas panel (e) shows a 
closeup. Full lines show the charge correlation function at 
time t after the quench, whereas dotted lines are results 
for free fermions. 
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FIG. 1. (Color online) Charge correlations S p (r, t) (full lines) 
at different times t after the quench. Here, A = 0.5, loo = 0.5J, 
L = 22 and /3J — 44. Dotted lines correspond to S p (r) for 
A = 0. Panel (e) shows a closeup of panel (d), and additional 
results for t = 0.63J. 



At t = 0, Fig. [TJa), we observe the clear 2/cp charge 
correlations expected for the Peierls state. At t > 0, af- 
ter the quench to A = 0, these oscillations decay over a 
time scale of about J / v2, after which 2fcp oscillations are 
practically absent at large distances. Shortly before that 
time, namely for t ss 0.63 J, small-amplitude 2kp oscilla- 
tions comparable to the free-fermion results are visible, 
see Fig. [lje). For t > 0.7 J, a revival towards a charge- 
ordered state begins. Whereas the charge correlations for 
t = 0.7J in Fig. [ljd) agree well with those for the non- 
interacting system on the scale of the figure, the closeup 
in Fig. [lje) reveals that noticeable differences remain at 
small distances even for the 'optimal time' of t = 0.63J. 

Figure [2] shows the corresponding results for the spin 
correlation function S a (r, t) . The times shown are the 
same as in Fig. fTj In contrast to S p (r,t), Fig.[lja), there 
are no discernible 2fcp correlations at t = 0, in accor- 
dance with the nonzero spin gap of the Peierls state. Af- 
ter the quench, the differences from the noninteracting 
results visible at short distances diminish as a function 
of time. For t = 0.7 J, similar to the charge correlations, 
S a (r,t) looks very similar to the noninteracting result on 
the scale of Fig. [2jd), but the closeup in Fig. [2^e) again 
reveals clear deviations at small distances. Interestingly, 
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FIG. 2. (Color online) Spin correlations S a (r, t) (full lines) at 
different times t after the quench. Here, A = 0.5, uuo = 0.5J, 
L — 22 and f3J = 44. Dotted lines correspond to S a (r) for 
A = 0. Panel (e) shows a closeup of panel (d), and additional 
results for t = 0.63J. 
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FIG. 3. (Color online) Double occupation as a function of 
time for A = 0.5 and two values of the phonon frequency. 
Here, L — 22 and f3J = 44. The vertical line illustrates the 
noninteracting value 0.25. 



icantly smaller double occupation at t — for cj = 5J 
than for ujq = 0.5J. At t > 0, we observe a fast decay of 
the double occupation, with periodic oscillations around 
the noninteracting value of 0.25. The time scale associ- 
ated with this decay agrees with the decay of 2kp charge 
correlations in Fig. [T] A partial revival to about half the 
value at t — is visible at time t 5J, and reappears 
with this period also at later times (data not shown). 



CONCLUSIONS 



the results for S a (r,t — 0.7 J) in Fig. u\e) agree within 
the symbol size used with those for Sp(r,t = 0.7 J), de- 
spite substantial differences at t = 0. 

Because the system is quenched to a noninteracting 
and hence integrablc Hamiltonian (Ho), real thermaliza- 
tion is not possible!^ Instead, the time evolution seen 
in Figs. [T] and [2] is related to dephasing. Indeed, we ob- 
serve the expected periodic but incomplete revivals of the 
system at longer times, see Fig. [3] and discussion below. 
Conserved quantities such as the total energy and the mo- 
mentum distribution function are independent of t. Since 
the evolution is determined by Ho [cf. Eq. Q], there is 
only a single energy scale, namely the hopping integral 
J. However, we can modify the initial state by varying 
the parameters A and ujq. As emphasized before, our 
method fully accounts for quantum phonon fluctuations 
in the initial state. An increase of the phonon frequency 
brings the initial state much closer to the Peierls phase 
boundarypO To illustrate the dependence of charge order 
and its time evolution on the phonon frequency, we show 
in Fig. [3] the double occupancy, (n^nn) , which reflects 
the pairing of electrons at the same site in the Peierls 
state, and its decay after the quench. Comparing the 
cases of oj — 0.5 J and ujo — 5 J, Fig. [3] reveals a signif- 



Motivated by an increasing number of experiments 
on ultrafast dynamics after photo-induced phase tran- 
sitions, and by the enormous challenge of numerically 
describing time-dependent phenomena in systems with 
electron-phonon interaction, we considered the quench of 
a Peierls insulator to a noninteracting Hamiltonian. Us- 
ing an extension of the continuous-time quantum Monte 
Carlo method, we presented exact results for equal-time 
charge and spin correlation functions. The time evolution 
of the latter is related to dephasing, and is determined by 
the only remaining energy scale after the quench, namely 
the hopping J. We find that the dominant 2/cf charge 
correlations of the initial Peierls state decay over a time 
scale of about J/v2- After that time, we find that spin 
and charge correlations are virtually identical. 

In the present case, the electron-phonon coupling and 
the quantum nature of phonons only play a role for the 
initial state. Increasing (decreasing) the phonon fre- 
quency (the coupling) moves the system closer to the 
metallic state, thereby suppressing charge order. Impor- 
tantly, the results shown cannot be obtained by consider- 
ing the antiadiabatic limit ujo ~ * o°j where the Holstein 
model maps onto an attractive Hubbard model. The rea- 
son is that the attractive Hubbard model at half filling 
always remains metallic, with the Luttinger liquid pa- 



5 



rameter K p fixed to one by symmetry.^ 

The exact results presented here can serve as bench- 
marks for more realistic calculations using, for example, 
the density-matrix renormalization group. In the latter, 
the phonon Hilbert space has to be truncated for the ini- 
tial state as well as for the time evolution. Finally, our 
findings may even become directly relevant for experi- 
ments with cold atoms or trapped ions; several different 
proposals exist how to r ealize el ectron-phonon Hamilto- 
nians in such systems!^ 135 1 37 1 38 1 
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